Melt-growth dynamics in CdTe crystals 
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We use a new, quantum-mechanics-based bond-order potential (BOP) to reveal melt-growth dy- 
namics and fine-scale defect formation mechanisms in CdTe crystals. Previous molecular dynamics 
simulations of semiconductors have shown qualitatively incorrect behavior due to the lack of an 
interatomic potential capable of predicting both crystalline growth and property trends of many 
transitional structures encountered during the melt — > crystal transformation. Here we demonstrate 
successful molecular dynamics simulations of melt-growth in CdTe using a BOP that significantly 
improves over other potentials on property trends of different phases. Our simulations result in 
a detailed understanding of defect formation during the melt-growth process. Equally important, 
we show that the new BOP enables defect formation mechanisms to be studied at a scale level 
comparable to empirical molecular dynamics simulation methods with a fidelity level approaching 
quantum-mechanical methods. 



CdTe-based materials have been instrumental for tech- 
nological breakthroughs in the semiconductor industry. 
The largest science and technology impact has been their 
widespread use in solar cells, radiation detectors, and 
medical imaging devices. While CdTe solar cells cur- 
rently have the lowest cost compared to any other photo- 
voltaic technologies pQ , the material is far from optimum 
as the record energy-conversion efficiencies achieved to- 
day are only 16%, significantly below the 29% theoretical 
value [2j[3]. CdTe-based Cdi^Zn^Te (CZT) alloys are 
currently the leading semiconductors for 7-ray detection, 
but their application is limited by low manufacturing 
yield (and, therefore, high cost) of detector-grade ma- 
terials. The under- achievement of CdTe solar cells and 
CZT detectors can both be attributed to charge-trapping 
defects formed during CdTe growth [IHS]. 

Direct molecular dynamics (MD) simulations provide 
a fundamental understanding of the CdTe growth dy- 
namics and defect formation. However, such simulations 
are extremely challenging because they sample a large 
number of metastable configurations not known a pri- 
ori. If the interatomic potential used in a simulation 
over-stabilizes a metastable configuration, that configu- 
ration will likely persist, leading to an amorphous growth. 
Previous MD simulations of semiconductor growth have 
shown qualitatively incorrect behavior due to the lack of 
an interatomic potential capable of predicting both crys- 
talline growth and property trends of many transitional 
structures encountered during the growth. In addition, 
past MD simulations of semiconductor crystalline growth 
were limited to vapor deposition [THE], wn ^ e cases for 
melt-growth (which are likely to be more challenging be- 
cause the crystallization occurs from a condensed liquid 
phase rather than the atom-by-atom assembly during va- 
por deposition) have yet to be demonstrated. 

A vast majority of successful MD simulations of semi- 
conductor vapor deposition [7HTT] used Stillinger- Weber 
(SW) [14 potentials. We established previously [15 that 
while SW potentials can easily ensure the lowest energy 



for the equilibrium tetrahedral semiconductor crystal and 
its crystalline growth during vapor deposition, they can- 
not satisfactorily capture the property trends of other 
configurations. Hence, they cannot accurately reveal de- 
fect formation. Tersoff potentials [16], on the other hand, 
capture property trends more accurately. However, there 
is no obvious way to ensure the lowest energy for the 
tetrahedral structure. Because the tetrahedral structure 
must have a lower energy than any other structures, it is 
unclear which phases should be included in the potential 
parameterization. As a result, crystalline growth is diffi- 
cult to achieve with Tersoff potentials unless the poten- 
tial parameterization is done iteratively with crystalline 
growth used as a criterion to actively select and modify 
the phases to be included in the fitting. Not surprisingly, 
we found that many literature Tersoff potentials [T7UT9] 
predict amorphous growth during vapor deposition sim- 
ulations. 

The objective of this work is to fill the missing re- 
search areas identified above by not only demonstrat- 
ing the crystalline growth of semiconductor compounds 
from melt, but also advancing beyond Tersoff potentials 
on capturing properties of different phases using a new 
bond-order potential (BOP) [2QH23] . This work will be- 
gin to allow detailed investigations of defect formation 
mechanisms during melt- or vapor- phase growth. While 
we focus on CdTe compounds, the methods are applica- 
ble to a broad range of semiconductors. 

Unlike the Stillinger- Weber [14] and Tersoff /Brenner 
[16] [24] potentials commonly used for semiconductors, 
BOP is analytically derived from quantum mechani- 
cal theories [20-23J. In particular, It incorporates pri- 
mary (<j) and secondary (tt) bonding and the valence- 
dependence of the heteroatom interactions, with the 
functional forms of the potential derived from tight- 
binding theory under the condition that the first two 
levels of the expanded Green's function are retained. De- 
tails of the parameterization of the CdTe BOP are dis- 
cussed elsewhere [25 . This parameterization considers 
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properties of a variety of elemental and compound con- 
figurations (with coordination from 1 to 12) including 
clusters, bulk lattices, defects, and surfaces, in addition 
to the crystalline growth of vapor deposition. 

There are currently two other CdTe interatomic po- 
tentials available in the literature, one [26 uses the 
Stillinger- Weber format [14], and the other [17] uses a 
Rockett modification [27 of the Tersoff format (TR). 
To evaluate different potentials, cohesive energies of a 
variety of Cd, Te, and CdTe phases calculated using 
various models are compared with the corresponding 
values obtained from our high-level density functional 
theory (DFT) calculations in Fig. [I] Here, various 
lattices are abbreviated as diamond-cubic (dc), simple- 
cubic (sc), body-centered-cubic (bcc), face-centered- 
cubic (fee), hexagonal-close-packed (hep), graphite (gra), 
graphene (grap), 7-Se (A8), zinc-blende (zb), wurtzite 
(wz), NaCl (Bl), and CsCl (B2). For clarity, these struc- 
tures are sorted to give monotonic DFT energy trends. In 
Fig. 1, the unfilled stars show the experimental cohesive 
energies [28] of the equilibrium phases, while the straight 
lines connecting the neighboring data points merely guide 
the eye. Because the DFT calculations typically give ac- 
curate energy trends but not the absolute energies, the 
cohesive energies obtained from DFT calculations are 
scaled to match the experimental values for the equilib- 
rium phases. Fig. 1 indicates that the cohesive energies 
calculated from BOP (solid lines) slightly oscillate (i.e., 
are non-monotonic) around the DFT benchmarks (thick 
light lines); however, these variations are quite minor. 
In fact, the BOP energy trends are considerably closer 
to those predicted by DFT than the corresponding re- 
sults of the SW and TR parameterizations. Most impor- 
tantly, BOP correctly specifies the lowest energies for the 
equilibrium phases of both elements and the compound, 
namely, the hep Cd, the A8 Te, and the zb CdTe, and the 
calculated cohesive energies of the lowest energy phases 
also match the corresponding experimental values. In 
sharp contrast, the lowest energy phases are calculated 
to be dc Cd, dc Te, and zb CdTe by the SW parame- 
terization and dc Cd, bcc Te, and B2 CdTe by the TR 
parameterization, with SW having the only correct result 
of zb CdTe. These results indicate that the TR parame- 
terization cannot be used to study any of the equilibrium 
Cd, Te, and CdTe phases as the structures will not even 
be stable in MD simulations. While the SW parameteri- 
zation can be used in some sort of MD simulation to study 
the equilibrium CdTe phase, caution should be taken in 
explaining the results involving defects as the potential is 
not transferrable to Cd and Te (and hence the defective) 
regimes. As a result, our new CdTe BOP approach sig- 
nificantly improves over other widely-used potentials on 
energy trends of different configurations leading to better 
description of defects. 

The geometry of the melt-growth MD simulations is 
shown in Fig. [2] An initial zb CdTe crystal containing 
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FIG. 1: Cohesive energies of a variety of Cd, Te, and CdTe 
phases calculated by various models. 



7800 Cd and 7800 Te atoms with L = 260 (400) layers 
(about 450 A) in the x- direction, H = 40 (044) layers 
(about 50 A) in the y- direction, and W = 24 (044) layers 
(about 30 A) in the z- direction was first created using 
the equilibrium lattice constant. Periodic boundary con- 
ditions were used in all three coordinate directions so 
that the system can be viewed as infinitely large. The 
two ends in the x- direction containing 2 (400) planes 
(about 1.8 A) were maintained at K temperature (i.e., 
atom positions were fixed) so that these two regions acted 
as seeding crystals. The two regions containing 24 (400) 
planes (about 41 A) adjacent to the fixed ends were con- 
trolled at a low temperature of T\ ow = 1000 K. A region 
containing 128 (400) planes (about 220 A) in the middle 
of the sample was controlled at a high temperature of 
Thigh- The remainder of the sample was left free. The 
middle portion of the sample was then melted by first 
giving a random displacement to all the atoms not in 
the fixed regions, and then annealing the system using a 
MD simulation with T high > 2200 K. The melt-growth of 
the crystal was then simulated in a second MD run at a 
different desired Thigh temperature where the Ti ow and 
Thigh regions grew/shrank at a nominal growth rate of R 
= 0.2 A/ps, as shown in Fig. [2] (the accelerated growth 
rate is required for MD simulations). 

Simulations were performed at two different high tem- 
peratures of Thigh = 2200 K and T high = 1800 K. The 
projected x-y configurations are shown in Fig. [3] as a 
function of time. The temperature profiles along the x- 
direction are superimposed on the atomic configurations. 
Fig. |3ja) shows the configuration prior to the growth 
where the middle section of the sample was melted at a 
temperature of Thigh — 2200 K. Fig. [3^b) shows that 
at t = 0.3 ns, the crystal/melt interface moved to a lo- 
cation corresponding to a temperature between 1300 
K and 1500 K in both the T high = 2200 K and T high = 
1800 K growth conditions, in good agreement with the 
melting temperature of CdTe, T m = 1365 K. In contrast 
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FIG. 2: Geometry of molecular dynamics simulations of melt- 
growth. 

(a) t = 0.0 ns 



— 1500 



50 
(b) t = 0.3 ns 

2500 




200 250 300 350 400 450 
x (A) light balls: Cd; dark balls: Te 

FIG. 3: System configurations (projected on the x-y plane) 
and temperature profiles (along the x- direction), (a) prior 
to growth; (b) 0.3 ns after growth started using two growth 
conditions of T high = 2200 K and T high = 1800 K; and (c) 
0.9 ns after growth started using Thigh = 2200 K and Thigh 
- 1800 K. 



to conventional SW or TR-based potentials, our BOP- 
based MD method allows a physically-correct crystalline 
growth from the melt while at the same time predicts 
accurate energy trends of various configurations. 

Fig. |3|c) indicates that the temperature profiles 
dropped to 1000 K across the sample length at t = 0.9 ns 
for both simulation conditions. Interestingly, the sample 
became entirely crystalline at Thigh — 1800 K but the 
middle portion remained amorphous at Thigh = 2200 K. 
This amorphous zone did not change with further simu- 
lation at 1000 K, but was found to crystallize when the 
temperature of the middle portion was raised to 1500 
K and then slowly cooled to 1000 K. Fig. [3J therefore, 
reveals how defects are trapped at a high growth rate. 



To further explore defects in the grown crystals, the 
crystallinity parameter developed previously [12] was cal- 
culated along the x- direction for the two configurations 
shown in Fig. [3^c), and the results are shown in Figs. 
|4ja) and|4jb) respectively. Fig. [4ja) shows that the crys- 
tallinity for the sample obtained from the Thigh — 

2200 K 

condition is uniform along the sample length except be- 
tween 170 A and 280 A where the crystallinity exhibits a 
drop, in agreement with the trapped amorphous phase in 
this region. Surprisingly, the sample obtained at Thigh — 
1800 K also exhibits a crystallinity drop near the center 
of the sample even though no amorphous zone is iden- 
tified. To understand this, the low crystallinity region, 
which is framed in Fig. |4^b), is magnified and examined 
in Fig. [5|a). Fig. [5|a) indicates a planar defect where the 
projected Cd — » Te stacking, shown as an arrow in the 
negative x- direction, is rotated by about 109.47° clock- 
wise. Such a rotation can originate from a [21 1] /6 slip 
on the (111) plane. The defect, therefore, corresponds to 
a stacking fault, and causes the drop in crystallinity ob- 
served in Fig. [4^b). Similar stacking faults have been ob- 
served in experiments [29j [30] . A time-resolved analysis 
surprisingly indicates that the growth front is not always 
perpendicular to the growth direction; rather, it forms 
local trailing {111} facets, suggesting a higher stability 
and a slower solidification on these planes. In particular, 
the formation of such facets was often accompanied by 
the nucleation of stacking faults, Fig. |5|b). This discov- 
ery is further supported by a previous MD study on gold 
melt-growth, where {111} planes were found to cause 
stacking faults and slow solidification kinetics [31] . While 
mechanical twinning is known to be caused by stresses, 
our simulations provide a mechanistic explanation for the 
formation of a stacking fault: when growth occurs on a 
hexagonal {111} plane, say, plane A, atoms can occupy 
two different sites B and C. One of these sites corresponds 
to the lattice sites and the other corresponds to stacking 
fault sites. The energy difference between these two sites 
is small, and hence, it is not surprising that the conden- 
sation on a {111} plane may nucleate a stacking fault 
defect. Since the (100) and (110) planes have only one 
low energy (lattice) site, such stacking faults do not form 
if growth occurs strictly on (100) or (110) planes. These 
results suggest that using {100} or {110} growth planes 
can reduce stacking faults if the growth technique per- 
mits a sufficiently high temperature gradient and a suffi- 
ciently low growth rate to prevent the formation of local 
{111} interfaces. However, if the growth technique does 
not permit a sufficiently high temperature gradient and 
a sufficiently low growth rate, the {111} growth planes 
might be beneficial despite its high propensity for stack- 
ing faults because this interface is more likely to remain 
fiat to minimize dendritic growth. 

In conclusion, we demonstrate that (a) new BOP-based 
MD simulations can accurately predict melt-growth of 
semiconductors; (b) strictly derived from a quantum- 
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FIG. 4: Crystallinity analysis of the configurations shown in 
Fig. [gc). (a) T h i gh = 2200 K and (b) T high = 1800 K. The 
framed region will be further analyzed in Fig. [5] 
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FIG. 5: (a) Defect analysis of the framed region shown in 
Fig. |2Jb); and (b) Growth front at the time when the same 
stacking fault shown in Fig. |5ja) was nucleated. Dashed line 
shows stacking fault, and thick line represents growth front. 



mechanical formalism, BOP enables defect formation 
mechanisms to be studied at a scale comparable to em- 
pirical MD methods and a fidelity approaching quantum- 
mechanical methods; and (c) amorphous defects can be 
trapped, and stacking faults can nucleate on {111} facets 
surprisingly formed during the melt-growth of semicon- 
ductors in non <111> directions. 
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